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Abstract 

We study a stochastic process describing the onset of spreading dynamics of an epidemic in a 
population composed by individuals of three classes: susceptible (S), infected (I), and recovered 
(R). The stochastic process is defined by local rules and involves the following cyclic process: 
S^I^R^S (SIRS). The open process S^I^R (SIR) is studied as a particular case of the SIRS 
process. The epidemic process is analyzed at different levels of description: by a stochastic lattice 
gas model and by a birth and death process. By means of Monte Carlo simulations and dynamical 
mean-field approximations we show that the SIRS stochastic lattice gas model exhibit a line of 
critical points separating two phases: an absorbing phase where the lattice is completely full of 
S individuals and an active phase where S, I and R individuals coexist, which may or may not 
present population cycles. The critical line, that corresponds to the onset of epidemic spreading, is 
shown to belong in the directed percolation universality class. By considering the birth and death 
process we analyse the role of noise in stabilizing the oscillations. 

PACS numbers: 87.10.Hk, 05.40.-a, 02.50.Ey 
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I. INTRODUCTION 



Among the several theoretical descriptions of population biology systems, the stochastic 
approach is one capable of revealing the ubiquitous fluctuations observed in these systems 
Using a stochastic lattice gas model, we study here the dynamics of propagation 
of an epidemic in a population in which the individuals are separated into three classes 
determined by their relative states to a given disease: susceptible (S), infected (I) and 
recovered (R) individuals. One of the most well known models in this context is the so- 
called susceptible-infected-recovered (SIR) model [10|, [ill, ll2|, [13|, 



a model for an epidemic which occurs during a time interval that is much smaller then 
the lifetime of the host. It describes the spreading of an epidemic process occurring in 
a population initially composed by susceptible individuals that become infected by contact 
with infected individuals. Once infected the individuals can recover spontaneously becoming 
immune. An extension of the SIR model, to be considered here, is obtained by allowing the 
recovered individual to becoming susceptible. This extension is called susceptible-infected- 
recovered-susceptible (SIRS) model and will be the main concern of this paper. The SIR 
model will be then a particular case of the SIRS model. 

The stochastic dynamics is studied on a two-dimensional lattice that represents the space 
where the individuals live. Each site of the lattice is occupied by just one individual that can 
be in one of the three states: susceptible (S), infected (I) or recovered (R). The stochastic 
rules are such that S becomes I by a catalytic reaction, I turns into R spontaneously, and R 
becomes S also spontaneously. By means of Monte Carlo simulations, dynamic mean-field 
approximations and also by means of Fokker-Planck and Langevin equations, we study the 
dynamics and the onset of an epidemic spreading. Our results show that the stationary 
states described by the SIRS stochastic lattice gas model, can be of two types: an absorbing 
(or inactive) state, where the lattice is fully covered by susceptible, and an active state. In 
the active state the processes of infection, recovering and loss of immunity are continuously 
occurring. The phase transitions that take place in this model are studied in detail, as well 
as the critical behavior, and critical exponents are obtained from Monte Carlo simulations. 
We remark that the transitions are nonequilibrium ones because the SIRS stochastic lattice 
gas model belongs to the class of models that are intrinsically irreversible. In other words, 
the dynamics do not obey detailed balance. 
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From the stochastic lattice gas model we derive a birth and death process for the numbers 
of individuals of each class and analyze relevant quantities to characterize the epidemic 
spreading such as the time correlation functions. We find that the SIRS model presents 
stochastic oscillations on a certain region of the space of parameters, inside the active region. 
These oscillations are better understood if the system is described by a Langevin equation 



from which it becomes clear that the amplitude of these oscil 



increases and that the noise enhances the oscillations jll. 14. 15. 16. 17|. The Fokker-Planck 
and the associate Langevin equations are derived from the birth and death process which 
in turn is obtained from the original SIRS stochastic lattice gas model. The derivation of 
the birth and death process is obtained by a contraction procedure in which the degrees of 



ations decay as the system size 



freedom is reduced to a few variables 



We have also analyzed the critical behavior of the model occurring around the transition 
from the active to the nonactive state. We find that the critical behavior places the SIRS 
but not the SIR model into the directed percolation (DirP) universality class. The critical 
behavior is obtained by determining the static and dynamic critical exponents by means of 
numerical simulations. The SIR stochastic lattice gas model studied here is shown to belong 
to a distinct class, namely the dynamic percolation (DynP) universality class in agreement 
with the conjecture by Grassberger [8(. 



II. STOCHASTIC LATTICE GAS MODEL 



The SIRS stochastic lattice gas model is defined on a regular lattice of N sites in which 
each site can be occupied by either one susceptible individual (state S), or one infected 
individual (state I), or one recovered individual (state R). The dynamics consists of three 
subprocesses: (1) I + S — > I + I, auto-catalytically; (2) I — > R, spontaneously; and (3) R — > 
S, spontaneously. At each time step a site is randomly chosen, (i) If the site is in the state 
S and there is at least one neighboring site in the state I then S becomes I with probability 
proportional to a parameter b and the number of neighboring I sites; (ii) if the site is in 
state I it becomes R spontaneously with probability c; and (iii) if the site is in state R it 
becomes S spontaneously with probability a. 

At each site % of a square lattice we attach a variable t]i that takes the values 0, 1 or 2, 
according to whether the site is in the state R, S, or I, respectively. The allowed transitions 
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of the states of a site are cyclic, that is, — > 1 — > 2 — > and the corresponding transition 
rate is given by 

MV) = \&(iH, 1) E <%> 2 ) + 2) + 0), (1) 
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where r\ = (r]i, . . . , rji, . . . , t^at) denotes the whole microscopic configuration, 6(x, y) denotes 
the Kronecker delta, the summation is over the nearest neighbours of site % and b, c and a 
are parameters related to the subprocesses (i), (ii) and (iii) described above. The transition 
rates of the SIR model is regarded as a particular case of the SIRS transition rate given 
by equation ([I]), obtained by setting a = 0, which means that the subprocess (3) described 
above is suppressed. 

The time evolution of the probability distribution P(f], t) at time t obey the master 
equation 

j t P(rj, t) = ^{wtMPtf, t) - w^Pirj, t)}, (2) 

where the configuration if is obtained from configuration rj by an anticyclic permutation of 
the state of site i, that is, 2 — > 1 — > — > 1. 

By rescaling time in the master equation it is always possible to reduce the three parame- 
ter of the model to just two. Thus we choose the parameters in such a way that a + b + c = 1 
so that a, b and c can be interpreted as probabilities; this allows us to introduce a new 
parameter p such that a and b are given by 

1 — c 1 — c , . 
a = — p and b = — h p, (3) 

where -1/2 < p < 1/2 and < c < 1. 

We have obtained the phase diagram of the SIRS stochastic lattice gas model defined on 
a regular square lattice by using Monte Carlo simulations. The phase diagram in the space 
of parameters p and c is shown in figure HJ The model exhibits an absorbing susceptible 
state and an active state where individuals are continuously changing their states. We 
have located the critical line separating the two phases and determined the dynamic and 
static critical exponents. Both are in agreement with the ones of the direct percolation 
universality class. The dynamic exponents were obtained by time dependent Monte Carlo 
simulations [18] and the static exponents by Monte Carlo simulations performed on finite 
square lattices with a number of sites up to 256x256. In the last case we have used a finite 



size scaling proposed in reference 



191 ] . Our best results for the exponents are given in table 
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FIG. 1: (Color online) Phase diagram of the SIRS stochastic lattice gas model in the p-c space of 
parameters showing the two phases: a nonactive phase, with only susceptible individuals, and an 
active phase, with the three types of individuals, where the disease spread over the population. The 
phase boundaries were obtained by Monte Carlo (MC) simulations, by pair mean-field (PMF) and 
simple mean-field (SMF) approximations. The right edge of the triangle, where a = 0, corresponds 
to the SIR model. 

HI Our estimates for the SIRS model (a > 0) are in agreement with the ones for the directed 



percolation universality class 



20). 



For the SIR stochastic lattice gas model (a = 0) we have used the same Monte Carlo 
procedures to determine the static and dynamic critical exponents. Our estimates of the 
exponents, shown in table HI are in agreement with the ones related to the dynamic percola- 
tion universality class |s[ 0]. We remark that the active phase of the SIR stochastic lattice 
gas model is characterized by infinitely many absorbing configurations. These absorbing 
configurations are reached when the number of infected individuals disappears. 

It is worth mentioning that when the values of the parameters b and c are much smaller 
than a and each other of the same order, corresponding to the region around the left corner of 
the triangular phase diagram of figure HJ an R individual is almost instantaneously converted 
into an S individual. We may therefore replace the two spontaneous processes I— >R and 
R— >S by just one spontaneous process I— >S. The SIRS model is thus reduced to the process 
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TABLE I: Critical exponents from Monte Carlo simulations for the SIRS model and SIR model 
(a = 0) denned on regular square lattice. For comparison we show the exponents related to the 
directed percolation (DirP) universality class and to the dynamic percolation (DynP) universality 



class according to reference 



20]. 



a 


c (3 (3/v r\ 


z 


5 


0.00 


0.1765(5) 0.129(10) 0.096(7) 0,586(1) 


1,773(1) 


0,094(3) 


DynP 


0.1389 0.1042 0.586 


1.771 


0.092 


0.025 


0.19640(5) 0.574(13) 0.806(11) 0.22(2) 


1.2(1) 


0.42(3) 


0.05 


0.21355(5) 0.56(2) 0.805(14) 0.232(7) 


1.163(7) 


0.424(15) 


0.10 


0.23270(5) 0.56(2) 0.83(3) 0.230(8) 


1.126(10) 


0.435(12) 


0.20 


0.2377(1) 0.56(2) 0.81(2) 0.227(2) 


1.132(2) 


0.446(2) 



DirP 0.583(4) 0.795(10) 0.2295(10) 1.1325(10) 0.4505(10) 



S— >1— ►S, composed 
models called SIS 



}y a catalytic infection and a spontaneous recovering, denning a class of 



211 ] ; in the present case it corresponds to the contact process, known 
to be in the DP universality class [18| . In the limit of small values of b and c, the ratio b/c 
is identified as creation rate A of contact process, so that the transition line around the left 
corner of the triangular phase diagram is given by b/c = X c where A c is the critical infection 
rate of the contact process defined on a regular square lattice. 

III. BIRTH AND DEATH PROCESS 

A distinct stochastic description of the system under discussion is the one known as birth 



and death process or one-step process |22j. In this description the state of the system is 
characterized by a few stochastic variables. A birth and death process is here obtained 
from the SIRS stochastic lattice gas model by reducing the number of degrees of freedom. 
We use a contraction procedure leading to a description with just two stochastic variables: 
the number n of susceptible individuals and the number m of infected individuals. These 
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FIG. 2: Transitions of the birth and death stochastic process in the space of susceptible and infected 
numbers for the SIRS model. The transition to the east represents a recovered individual becoming 
susceptible with rate A + q, to the northwest a susceptible individual becoming infected with rate 

A |_, to the southwest an infected individual becoming recovered with rate Ao_. The full circle 

represents an absorbing state. 

variables are defined by 

N N 



n 



$3%, 1) and m = £%,2). (4) 



The appropriate stochastic description involves the probability P(n, m, t) of n and m at time 
t which is assumed to obey the following birth and death master equation 

d 1 1 

—P(n, m,t) — N {A aT {n — a,m — r)P(n — a,m — r,t) 

dt <j=-1t=-i 

— A aT (n,m)P(n,m,t)}. (5) 

This equation comes from the master equation ([2]) by summing over 77 with the restriction 
that the number of susceptible individuals is n and the number of infected individuals is m. 
This procedure leads to the following results 

A +0 = az, (6) 

A_+ = bxy, (7) 
A)- = cy, (8) 
7 



N 



n (susceptible number) N 

FIG. 3: Transitions of the birth and death stochastic process in the space of susceptible and infected 
numbers for the SIR model. The transition to the northwest represents a susceptible individual 

becoming infected with rate A |_ , to the southwest an infected individual becoming recovered with 

rate Ao_. A recovered individual never becomes susceptible. The full circles represent absorbing 
states. 

where z — 1 — x — y, x = n/N and y = m/N, and N is the population size. Expressions (J6]) 
and ([8]) are exact and expression ([7]) is derived by a scheme similar to the so called simple 



mean-field approximation 17]. 



The birth and death process defined by the master equation ([5]) can be regarded as 
a random walk in the space (n, m) as shown in Figure EJ The possible jumps are (a) 
(n, m) — > (n + l,m), with probability A +0 , (b) (n, m) — > (n — 1, m + 1), with probability 
A_ + , and (c) (n,m) — > (n,m — 1) with probability A _. 

We have simulated the birth and death stochastic process according to the transition 
rates given by equations ([6]), ([7]) and (El) for iV = 5000. We have also determined the time 
correlation functions as a function of the time lag t defined by 

c f9 (t)= f[f(t' + t)-f}[g(t')-g}dt', (9) 



where / and g can be any of the variables n or m. As an example, we show in figure H] the 
auto-correlation function for infected individuals and the cross-correlation between infected 
and removed individuals for the following values of the parameters: p = 0.25 and c = 0.40. 
The damped oscillations of the time correlation function indicates that the populations of 
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FIG. 4: (Color online) Infected-infected (inf-inf) time autocorrelation function and infected- 
removed (inf-rem) time cross-correlation function versus the time lag for p = 0.25 and c = 0.40. 
Results from numerical simulation of the birth and death process. 

individuals in each class oscillate in time with period of T = 72(3) Monte Carlo steps. This 
period was confirmed by Fourier analysis of the time correlation functions. The active phase 
can thus be separated into two regions as can be seen in figure [5j One region corresponding 
to sustained oscillations in population of individuals and the other without oscillations. 



IV. LANGEVIN EQUATIONS 



Considering the situation where the population size N is large we perform an expansion 
of the master equation (jHJ) in powers of 1/N 22|, |23J]. Up to second order in 1/N we get a 
Fokker-Planck equation for the probability density V(x,y,t), namely 

d d d 



where 



and 



/i(ar, y)=az- bxy, f 2 {x, y) = bxy - cy, 

D n (x,y) = az + bxy, 



(10) 

(11) 
(12) 
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FIG. 5: (Color online) Phase diagram of the SIRS model according to the birth and death process. 
Inside the active phase, below the nonactive phase, there are two regions corresponding to the 
oscillatory (osc) and to the nonoscillatory (non-osc) regimes. The line of separation between these 
two regimes is given by equation (|24p . 

Dn(x,y) = -bxy, (13) 
-D22 (x, y) =bxy + cy. ( 14) 

The Langevin equations associated to the Fokker-Planck equation (jTUI) are 

§=/,(*,„)+ 4-c.w. 



(15) 



£ = A(»,») + -4««>, 



(16) 

where d(t), with i — 1,2, are white Gaussian noise functions with zero mean obeying the 
properties 

WKjW)) = D tJ (x,y)6(t-t'). (17) 
Similar expansion procedure to study population biology systems has been used by other 



authors 
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V. MEAN-FIELD APPROXIMATION 



In the limit N — > oo the noise disappears and we are left with a deterministic dynamics 
which is identified here as a dynamic simple mean-field approach. The time evolution of the 
densities is given by 

= az - bxy, (18) 

% = hx v - °y- ( 19 ) 

The equation for z is not necessary because x + y + z = 1. The trivial solution is 
(x*,y*) = (1,0) corresponding to a state where the entire population is composed by sus- 
ceptible individuals only, called the inactive phase. The nontrivial solution is given by 

x* = c/b, y* = a(b-c)/b(a + c), (20) 

and is valid as long as b/c > 1, corresponding to a state of epidemic spreading, where the 
individuals of the three classes coexist. The stability of the solutions are obtained from the 
eigenvalues of the Jacobian matrix 

(-a - by* -a - bx* \ 
. (21) 
by* -c + bx* I 

For the trivial solution the eigenvalues are —a and b — c. Therefore this solution is stable 
as long as b/c < 1. The transition from the inactive phase to the active phase occurs along 
the line defined by6 = corc = (1 — 2p)/3, as shown in figure [1] and [51 For the nontrivial 
solution, corresponding to the active state, the eigenvalues A are given by 

H(-t^4 (22) 

where 

Therefore, when A < 0, the eigenvalues are complex and the solution approaches an asymp- 
totic stable focus. When A > 0, they are real and the solution approaches an asymptotic 
stable node. The transition from a focus to a node occurs when A = 0, that is, when 

a 2 (a + b) 2 = 4a(b-c)(a + c) 2 , (24) 

which defines the separation line between these two behaviors, as shown in figure [51 
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Let us find some particular solutions of the separation line ([21]). When a = b (or p = 0) 
this equation reduces to c 2 +4c— 1 = which gives c = \fh — 2 = 0.236068. When b = 0+1/2 
(or p = 1/4), equation (12411 gives two solutions: one is c = 1/2 and the other is the solution of 
the root of c 2 + (5/2)c - (1/8) = which gives c = (3^ - 5)/4 = 0.049038. Improvements 
at qualitative level in the phase diagram can be achieved by considering dynamic mean- 
field approximations which takes into account correlations between neighbors sites. This 



procedure has been largely used in the context of stochastic lattice gas models, also ca 
interacting particle systems, describing population biology systems [10l 
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24 



25 



led 



27 



26j, 



281 ] . By using a pair mean-field approximation (PMF) 25] we find for the present model 



a transition line between the active phase and the inactive phase which is shown in figure [TJ 



VI. CONCLUSIONS 



In this article we have examined a spatial-structured model describing the dynamic of a 
SIRS epidemic process. This model shows a continuous phase transition between a phase 
where the disease spread and other where it dies out in a short time. The phase diagram of 
the SIRS stochastic lattice gas model was obtained by Monte Carlo simulations and dynamic 
mean-field approximations. Our Monte Carlo results indicate that the SIRS stochastic 
lattice gas model exhibits a line of critical points that belongs in the directed percolation 
universality class. As far as we know this is the first determination of the universality class 
of the SIRS stochastic lattice gas model. 

We have also derived a birth and death process for the numbers of individuals in each 
population class. This process was derived from the stochastic lattice gas model in the 
scope of simple mean-field approximation. A Fokker-Planck equation as well as the associate 
Langevin equations were then derived from the birth and death process. The simulation of 
this process, which can be seem as a random walk in the space of the number of susceptible 
and infected individuals, shows that the active state is subdivided into two regions: one 
in which the numbers of individuals in each class oscillate in time and another where they 
remain constant in time. 

In summary, we have devised and analyzed the SIRS epidemic process at different levels 
of description: (a) the spatial, stochastic and local description, given by the SIRS stochastic 
lattice gas model; (b) the stochastic but zero dimensional birth and death process; and (c) 



12 



the deterministic level in which the space is not considered, given by the simple mean-field 
approximation where all individuals are well mixed. The fundamental level is the one based 
on stochastic lattice gas models also called interacting particle systems j^ . These spatial 
structured models are relevant when the diseases spreads over a region with few individuals 



ated, loosing contact with the 



or when a group of individual of a given class becomes iso 
other class of individuals during a certain period of time [29 
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